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ABSTRACT 


Development  and  Validation  of  a  Vertically 
Two-Dimensional  Mesoscale  Numerical  Model.  (August  1985) 
Michael  Kent  Walters,  B.S.,  Texas  A&M  University 
Chairman  of  Advisory  Committee:  Dr.  Dusan  Djuric 


A  vertical,  two  dimensional,  grid-point  mesoscale 
model  is  developed,  using  the  equations  of  motion  and 
thermodynamics  in  a  dry  flow.  A  non-dimensional  vertical 
coordinate  s  is  used.  The  hydrostatic  assumption  is  made. 
To  avoid  the  sensitivity  of  the  continuity  equation, 
several  derived  equations  are  used  based  on  the  first  law 
of  thermodynamics,  the  pressure  tendency  equation,  and  the 
top  boundary  condition,  which  was  that  the  vertical  motion 
was  zero  at  the  top  boundary.  These  derived  equations 
include  the  equation  for | prediction  of  density, 
Richardson's  equation  for  vertical  motion,  and  the 
pressure  tendency  at  the  top  of  the  model.  A  simplified 
calculation  of  dry,  subgrid  convection  is  made  to  prevent 
instability.  The  equations  are  integrated  using  explicit 
finite  differences  with  a  time  step  of  2  minutes.  ^ 

A  test  made  with  a  static  start  and  no  forcing 
showed  no  change  in  initial  variables.  The  kinetic  energy 
budget  for  the  model  was  calculated  by  modifying  the 


initial  state  through  differential  heating.  The  test 


.‘-V 


showed  a  diffusion  of  kinetic  energy  of  less  than  1%  per 
hour.  A  calculation  of  mass  continuity  using  solid 
lateral  boundaries  showed  mass  continuity  preserved  to 
computer  accuracy.  Problems  which  occurred  in  duplicating 
the  low-level  jet  calculations  of  Shen  (1980)  were  solved 
by  rederiving  Richardson's  equation  with  a  top  boundary 
condition  in  which  the  vertical  gradient  of  vertical 
motion  vanished  at  the  top  of  the  model,  and  by  specifying 
a  temperature  gradient  at  the  left  boundary.  Results  were 
consistent  with  those  reported  by  Shen.  The  conclusion  is 
made  that  the  model  is  formulated  correctly  and  is  capable 
of  physically  realistic  results  in  two  dimensions  on  a 
vertical  plane..  These  results  are  consistent  with  the 
imposed  boundary  conditions. 
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1. 


INTRODUCTION 


Since  the  advent  of  computers  in  the  1940's  it  has 
become  increasingly  popular  to  apply  them  to  the  numerical 
solution  of  complicated  fluid  dynamics  problems  for  which 
no  analytical  solutions  are  possible.  One  of  the  original 
applications  of  computers  to  the  solution  of  such  problems 
was  in  the  field  of  meteorology.  Numerical  models  have 
been  developed  to  simulate  many  atmospheric  processes, 
from  land-sea  breezes  to  global  general  circulation.  With 
computer  resources  becoming  increasingly  available  and 
powerful,  many  investigators  have  become  interested  in 
applying  them  to  the  simulation  of  meteorological  events 
on  horizontal  scales  of  less  than  1000  km,  i.e.  mesoscale 
events . 

It  has  been  recognized  by  many  recent  investigators 
that  the  study  of  the  mesoscale  environment  is  very 
important  in  understanding  the  growth  and  development  of 
convective  clouds.  Ulanski  and  Garstang  (1978)  have  noted 
that  the  size  of  the  area  of  surface  convergence  is  a 
critical  factor  in  attempting  to  forecast  the  total  amount 
of  rainfall  which  can  be  expected  to  be  produced  by  a 
given  storm.  Tripoli  and  Cotton  (1980)  pointed  out  that 
the  level  of  kinetic  energy  which  an  individual  cumulus 
cloud  can  achieve  is  a  function  of  the  pre-existing 

The  citations  on  the  following  pages  follow  the  style 
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low-level  mesoscale  convergence  of  moist  static  energy. 

It  has  become  increasingly  evident  that  events  on  the  meso 
scale,  or  horizontal  scales  of  20  km  to  200  km  (Orlanski, 
1975) ,  are  of  primary  importance  in  the  study  of 
convective  phenomena.  For  this  reason,  it  is  desirable  to 
search  for  suitable  research  tools  to  study  these 
phenomena. 

This  study  was  undertaken  to  aid  in  the  development 
of  a  suitable  numerical  tool  which  could  be  used  to  study 
such  important  mesoscale  phenomena. 


2.  BACKGROUND  AND  STATEMENT  OF  THE  PROBLEM 

a.  Previous  studies 

As  part  of  the  Texas  High  Plains  Experiment,  or 
HIPLEX,  Djuric  and  Das  (1982)  recognized  the  need  for  a 
simple  operational  model  for  predicting  mesoscale 
convergence  fields  of  mass,  moist-static  energy  and 
vorticity.  They  realized  that  such  a  model  would  be 
limited  in  scope  but  expected  that  it  would  prove  powerful 
and  useful  in  the  pre-storm  and  early-storm  conditions 
when  the  feedback  from  the  clouds  to  the  environment  is 
relatively  unimportant.  The  existing  models  such  as  the 
NCAR-Drexel  model  (Perkey,  1976),  the  Florida  sea  breeze 
model  (Pielke,  1974;  Pielke,  1976)  ,  the  mesoscale  models 
developed  at  Pennsylvania  State  University  (Anthes  and 
Warner,  1978;  Warner  et.  al . ,  1978)  and  the  NOAA/ERL 
models  (Nickerson  et.  aK  ,  1978;  Nickerson  et.  <al.  ,  1979; 
Fritsch  and  Chappell,  1980)  were  deemed  unsuitable  for 
such  an  application  because  they  were  complex  and  demanded 
computer  time  which  was  considered  excessive  for  a  simple, 
real-time  forecast  model  of  the  type  envisioned.  They 
also  rejected  the  models  of  Tapp  and  White  (1976),  Klemp 
and  Wilhemson  (1978),  Schlesinger  (1978),  Cotton  and 
Tripoli  (1978),  and  Tripoli  and  Cotton  (1980),  because 
they  were  even  more  complex  than  the  previously  mentioned 
models  and  were  thus  unsuitable  for  their  purposes.  Faced 


A  three-dimensional  numerical  model  designed  to  meet 
the  requirements  mentioned  above  was  developed  by  Djuric 
and  Das  (1982).  This  model  was  named  the  TAMU  Mesoscale 
Numerical  Model.  When  they  entered  the  validation  phase 
of  the  model,  problems  began  to  become  evident  with  the 
program  itself.  These  difficulties  were  manifested  as 
numerical  instabilities  which  occurred  unpredictably  in 
certain  grid  points  of  their  model.  They  traced  the 
origin  of  the  numerical  instabilites  to  certain  advective 
terms  in  the  model,  but  were  unable  to  resolve  the 
difficulties.  Because  of  these  persistent  numerical 
instabilities,  they  postponed  planned  modifications  to 
their  boundary  conditions  which  they  hoped  would  improve 
the  model,  and  focused  their  research  efforts  on  the 
removal  of  the  numerical  instabilities.  However,  the 
problem  has  remained  unsolved. 

The  objective  of  this  study  is  to  design  and  verify  a 
vertically  two-dimensional  mesoscale  numerical  model  based 
on  the  equations  of  motion,  continuity,  and  thermo¬ 
dynamics,  which  is  capable  of  simulating  profiles  of  wind 
and  thermodynamic  variables.  This  model  will  be  similar 


to  the  model  of  Djuric  and  Das  (1982) ,  except  that  it  will 
be  only  vertically  two-dimensional.  Numerical  testing  of 
this  vertically  two-dimensional  model  may  give  insight 
into  the  nature  of  the  numerical  instabilities  which 
plagued  the  TAMU  model  because  the  predictive  and 
diagnostic  equations  and  the  finite-difference  forms  used 
in  the  two  models  will  be  very  similar.  In  addition,  it 
will  be  possible  to  test  various  boundary  conditions  which 
could  be  later  used  in  a  three-dimensional  version. 
Successful  completion  of  the  simpler  model  should  make 
possible  the  further  systematic  development  of  a 
three-dimensional  version  which  can  fulfill  the 
requirements  of  Djuric  and  Das  (1982)  mentioned  above. 


3.  THE  NUMERICAL  MODEL 

a.  The  coordinate  system 

The  model  uses  a  terrain-following  vertical  co¬ 
ordinate  system.  The  vertical  coordinate  is  a  non- 
dimensional  height  s  with  s  =  0  at  the  ground  and  s  =  1 
the  top  of  the  model,  as  shown  in  Fig.  1.  The  relation 
ship  between  the  vertical  coordinate  s  and  the  height  z 
above  sea  level  is  given  by 


S  =  (  Z  -  ZQ)  /  (  Zt  -  Z0  )  =  (  Z  -  Z0)  /  H 

where  Zg  is  the  terrain  elevation  above  sea  level, 
is  the  height  above  sea  level  of  the  top  of  the  model, 
H  is  the  model  thickness  in  meters.  The  height  Z  at  a 
coordinate  surface  s  =  constant  is  a  space  dependent 
variable,  but  is  independent  of  time. 


Fig.  1.  Diagram  of  the  vertical  coordinate  system 


The  model  uses  the  basic  equation  set  used  by  Djuric 


and  Das  (1982)  with  the  exception  that  the  equation  for 
continuity  of  water  vapor  is  not  used  because  the  model  is 
dry.  The  equations  are  as  follows: 

dV  *  _»  l  Id  6V 

dt  =  -  f  k  X  V  -  £  VP  -  *  Vz  +  ^  Is  (PKm  '  <2) 

If  "  -  P9H  ,  (3) 

§£  +  V  .  pv  +  |s  (ps)  +  g  V  •  VH  =  0  ,  (4) 


dT 

dt 


(5) 


p  =  pRT 


(6) 


The  symbols  to  be  used  here  are  defined  as  follows: 

^  =  horizontal  wind  vector 

k  =  vertical  unit  vector 

s  =  vertical  component  of  wind 

x  =  horizontal  coordinate 

s  =  vertical  coordinate 

t  =  time  coordinate 

p  =  pressure 

p  =  density 

f  =  Coriolis  parameter 


g  =  acceleration  due  to  gravity 
z  =  height  above  sea  level 
H  =  thickness  of  the  model 
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K  =  eddy  diffusivity  of  momentum  =  5  m  s 
m  -1 

T  =  eddy  stress 

Cd  =  surface  drag  coefficient  =  0.15 
T  =  temperature 
Q  =  heating  rate 

0^=  specific  heat  at  constant  pressure 
Cv=  specific  heat  at  constant  volume 

t  -  ycv 

R  =  gas  constant  for  dry  air 

Equation  (2)  is  the  horizontal  equation  of  motion, 

(3)  is  the  hydrostatic  equation,  (4)  is  the  continuity 
equation,  (5)  is  the  first  law  of  thermodynamics,  and  (6) 
is  the  equation  of  state.  Equations  (2) -(6)  form  a 
complete  set  for  dry  flow  when  Q  and  Km  are  specified. 

To  solve  this  equation  set  numerically,  several 
derived  equations  are  used  (Djuric  and  Das,  1982).  These 
are  based  on  the  pressure  tendency  equation 

St  “  <5 t>l  -  v  '  (HP9>  ds'  ‘  PHPi^  1 

s 

+  gHps  ,  (7 

the  first  law  of  thermodynamics,  and  the  boundary 


condition  s  =  0  at  the  top  of  the  model.  The  derived 
equations  used  are  the  following: 


H--v.7P-*§e.-e  - 


,  D  -  V  •  Vd 


r  at'  1 


■ 

V  •  Htf 

ds 

Equation  (8),  used  for  the  prediction  of  density,  i 
derived  from  the  first  law  of  thermodynamics  in  order  to 
avoid  the  sensitivity  of  the  continuity  equation,  (9)  is 
the  Richardson  equation  for  vertical  motion,  and  (10)  is 
the  equation  for  pressure  tendency  at  the  top  of  the 
model.  The  function  D  used  above  is  the  integrated 
divergence  above  the  level  s  and  is  given  by 


D  =  g  V  • (pHtf )  ds ’ 
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A  x  is  equal  to  45  km,  while  A  s  is  equivalent  to  250  m  in 
the  simulations  to  be  presented  here.  The  horizontal  wind 
component,  u,  is  defined  on  both  the  left  and  right 
lateral  boundaries,  while  pressure  is  defined  on  both  the 
upper  and  lower  boundaries.  As  shown  in  Fig.  2.,  the 
variables  are  subdivided  into  triangular  regions  carrying 
the  same  computational  indicies,  to  simplify  the  notation. 

d.  The  lateral  boundary  conditions 

A  simplified  lateral  boundary  condition  of  u  =  0  is 
used  for  purposes  of  model  validation.  This  boundary 
condition  was  chosen  because  it  allows  the  finite 
difference  formulas  used  in  the  interior  of  the  model  to 
also  be  applied  at  the  lateral  boundaries.  Due  to  the 
centered  nature  of  the  horizontal  space  derivatives,  use 
of  the  interior  finite  difference  formulas  at  the  grid 
points  immediately  adjacent  to  the  lateral  boundaries 
requires  the  knowledge  of  values  of  pressure  and  density 
outside  the  model.  Use  of  this  boundary  condition  makes 
it  possible  to  set  the  pressure  and  density  immediately 
adjacent  to  the  boundary  on  the  outside  equal  to  that 
immediately  adjacent  to  the  boundary  on  the  inside.  This 
excludes  the  need  for  introducing  separate  forward  or 
backward  space  differencing  at  the  boundaries  which  would 
complicate  the  validation  of  the  model. 


Method  of  solution  of  tne  equations 


Solution  of  the  equations  is  accomplished  in  the 
following  manner.  The  horizontal  divergence  is  calculated 
from  known  wind  components  and  densities  for  each  level  of 
the  model.  Equation  (11)  is  then  solved  at  each 
appropriate  grid  point.  Next,  the  following  individual 
parts  of  (10)  are  evaluated. 

D  -  V  •  Vp 


(V  *  HV) 


Once  these  are  calculated,  (14),  (15),  and  (16)  are 

integrated  over  the  depth  of  the  model.  The  pressure 
tendency  at  the  top  of  the  model  is  then  found  using  (10) , 
following  which  the  Richardson  equation,  (9),  is  solved. 

At  this  point  the  density  tendency,  (3),  is  computed. 

Next,  the  nonlinear  components  of  the  equations  of 
horizontal  momentum  are  evaluated.  Once  these  values  are 
known  throughout,  the  pressure  at  the  new  time  level  is 
obtained  from  integration  of  (3)  using  the  pressure 
tendency  at  the  top  of  the  model  and  the  density  tendency 
previously  calculated.  Temperature  is  then  obtained 
diagnostically  from  (6)  using  the  new  values  of  pressure 


and  density.  The  temperature  profile  is  checked  for 
possible  dry-convective  adjustment,  following  which  the 
pressure  and  density  are  recalculated  hydrostatically  if 
necessary.  Finally,  the  u  and  v  components  of  the  wind 
are  evaluated.  The  computation  then  proceeds  to  the  next 
time  step  by  repeating  the  entire  numerical  process. 

f.  The  finite-difference  equations 

The  preceding  computational  process  is  accomplished 

using  the  following  explicit  finite-difference  equations, 

in  which  DX  refers  to  the  horizontal  grid  increment,  DS 

refers  to  the  vertical  grid  increment,  variables  i  and  k 

refer  to  spatial  indicies  as  shown  in  Fig.  2,  and  DT  is 

the  time  step.  The  time  step  used  was  120  s,  which  was 

chosen  to  satisfy  the  Courant-Friedrichs-Levy  (CFL) 

condition  for  computational  stability.  To  satisfy  the  CFL 

DT 

condition  DT  must  satisfy  c^  <  1,  where  c  is  tne 
wave  speed  of  the  fastest  wave  expected  in  the  model, 
which  in  this  case  is  a  sound  wave  with  c  equal  to  340  m 
s"1,  which  was  based  on  the  warmest  temperature  expected 
in  the  model.  Superscript  variables  n  and  n+1  refer  to 
the  present  and  future  values  respectively.  Other 
superscripts  refer  to  exponentiation  in  the  usual  manner. 
Subscript  1  indicates  the  variable  is  defined  at  the  top 
of  the  model.  The  remaining  variables  have  been 


previously  defined. 


The  finite  difference  equations  are  as  follows: 


V  *  H  p  V  (i,k)  =  H (i  +  1)  *  (  p  (i+l,k)  +  p(i,k)  )  * 

u (i+1 ,k>  -  H ( i )  *  (  p ( i , k)  +  p(i-l,k)  )  *  u(i,k)  / 


(1 


(  2  *  DX  )  =  DV(i,k) 

1  n=12 

(  v  .  H  p  V  )  ds'  ( i ,  k )  =  ^  DS  *  G  *  DV  (i,n) 

)  n=k 


=  DIV  (i,k)  ,  (1 

V-Vp  (i,k)  =  V.pV  (i  ,k)  -  p  V  •  V  (i  ,k)  =  (  (p  (i,k) 

+  p(i,k+l)  +  p(i+l,k+l)  +  p(i+l,k))  *  u(i+l,k) 

-  (  p(i,k)  +  p(i,k+l)  +  p ( i-1 , k+1)  +  p(i-l,k)  }  * 

u  ( i  / k )  )  /  (2  DX)  -  (p  ( i  ,k)  +  p  (i  ,k  +  l)  )  * 
(u(i+l,k)  -  u(i,kl) )  /  DX  =  VDP  (i,k)  ,  (2 


DIV  -  V  p  /  y  p  (i  ,k)  =  (  DIV  (i  ,k)  +  DIV  (i,k  +  l)  - 

VDP  (i,k)  )  /  (  1.4  *  (p(i,k)  +  p (if k+1) )  )  = 

DVP ( i , k )  ,  (2 

pg-Ufk)  =  Q  (i,k)  /  (C  T(i,k)  )  =  QCT  ( i ,  k )  , 

p  (2 


—  V  •  V  H  (i,k)  =  (  H  ( i  + 1 )  *  u ( i  +  1 ,k)  - 

H 

H  ( i )  *  u  ( i , k)  )  /  HP  (i)  =  DHV(i,k)  , 


I 


s 

(  (D 


V-VP  /  TP)  +  QCT  -  DHV) ds 1  )  (i,k)  = 


DS  *  (  DVP ( i , n-1 )  +  QCT ( i , n-1 )  - 
DHV ( i , n-1 )  =  DIS(i,k) 


(2 


(2 


I  ds/P  (i,k)  =  DSP { i ,  k)  =  /  2DS  /  (p(i,n)  +  p ( i , n— 1 ) )  , 
0  ^ °  (25 

(|P>1  =  1.4  *  DIS (i ,kel)  /  DSP  (i,kel)  , 


s  (i ,k)  =  SD  (i , k)  =  DIS ( i , k)  -  DIS(i,kel)  *  DSP(i,k)  / 


DSP  (i,kel) 


(|£)  ( i ,  k)  =  -p(i,k)  *  (  -DIS (i ,kel) / {  DSP(i,kel)  * 

(p (i,k+l)  +  p (i/k)  )  *  0.5  )  +  ( DVP ( i , k )  +  QCT(i,k)) 
- ( (SD (i ,k+l)  *  (p ( i , k+1)  -  p(i,k) )  +SD(i,k)  * 
(p(i,k)  -  p(i,k-l))  /  2  DS  +  (u(i+l,k)  *  (p(i+l,k)  - 
p(i,k))  +  u (i ,k)  *  (p(i,k)  -  p(i-l,k))  /  (  2  *  DX  ) 

(2( 

The  new  u  and  v  components  of  the  wind  are  computed 
using  a  finite  difference  scheme  which  uses  pressure  at 
the  new  time  level  and  trapezoidal  Coriolis  terms,  which 
are  averaged  between  the  old  and  the  new  time  levels. 
Although  the  use  of  u  and  v  at  the  new  time  level 
seemingly  makes  the  scheme  implicit,  it  is  possible  to 
solve  for  u  and  v  explicitly  using  the  following  set  of 
eguations . 


h*1  =  (  AU  *  a 

*  AV  ) 

/[p"*1  * 

(  1  *  a2  )]  ; 

(25 

W  =  (  AV  -  a 

*  AU  ) 

/[p"*1  * 

(  1  .  a2  )]  . 

( 3  ( 

In  these  expressions  A(J,  AV,  and  a  are  given  by  the 
following: 


AU  =  (  -gg  )“TX  *  DT  +  p“  *  (u“  +  a  *(  v“(i)  + 


)/2)  -  -|-x(pu2)n  -§-s(pus)n 


In»  * 


,pu  oh  ,n, 

(fr-ss  }  > 


DT 


(3 


AV  =  <  "  (puv)n  -  |-s(pvs)n  -  (puv/H  )n  ) 


*  DT  +  pn  *  (  vn  -  a  *  (  un  (i)  +  un (i+1)  )  /  2  ) 


-  <  i  fsPKm  S  >  *  DT  ' 


(3 


and 


a  =  f  *  DT  /  2 


(3 


The  internal  friction  terms  in  (31)  and  (32)  are 
evaluated  by  the  following  finite  difference  formulas  at 
all  levels  above  the  first  level. 

F  4  P*.  55  =  1  Kmli'k)  *  {  -  u,i'k>  1  - 

Km(i,k)  *  (  u(i,k)  -  u(i,k-l)  )  /  (DS2  *  H2)  ,  (3 

If  §5  PK,n  li  ’  <  Km(i'k)  *  1  v'i'k‘l>  -  ''(!■*>  )  - 

K  (i,k)  *  (  v ( i , k)  -  v(i,k-l)  )  /  (DS2  *  H2)  *  (3 

m 

Boundary  shear  stress  at  the  lowest  level  is 
calculated  from  the  following: 

^  *  Km(i,k)  *  (  u ( i , 2 )  -  u (i , 1)  )  /  DS2  -  Cd  * 

u(i,l)  *  u ( i , 1)  /  DS  (3 

H  li  =  Km(i,k)  *  (  v(i'1)  ~  vfi'1)  )  /  DS2  -  Cd  * 


v ( i , 1)  *  v ( i  ,  1 )  /  DS 


(3 


Eqs.  (36)  and  (37)  assume  that  the  eddy  stress  nT  is 
given  by  -j!jpKm  ^  and -^pKm  ^  in  the  layer  immediately 
above  the  lowest  layer,  while  the  eddy  stress  for  the 
lowest  layer  is  given  by  pCd  *  u  *|u|and  pCd  *  v  *  |v|  . 

is  multiplied  by  the  initial  density  at  the  beginning 
of  the  simulation  rather  than  the  updated  density  at  each 
future  time  step.  Similarly,  the  Coriolis  parameter  and 
the  surface  drag  coefficient,  which  are  assumed  constant, 
are  multiplied  by  the  initial  value  of  density  at  each 
grid  point  prior  to  the  first  forecast  step.  The  initial 
density  is  used  because  the  change  in  density  during  the 
simulation  is  small,  and  f,  Km  and  Cd  are  arbitrary 
constants.  Use  of  the  initial  density  for  these  . 
operations  takes  less  computation  time,  and  introduces 
small  error  since  the  constants  are  arbitrary. 

The  nonlinear  terms  in  the  above  formulas  are 
evaluated  according  to  the  following  finite  difference 
formulas : 

^  (pu2)  (i,k)  =  (  p(i,k)  *  (  u ( i , k)  +  u(i+l,k)  )2  - 
p (i-1 ,k)  *  (  u(i,k)  +  u ( i-1 ,k)  ) 2  )  /  (  4  *  DX  )  ; 

6  (38) 

g^(pus)(i,k)  =  (  (  p(i,k)  +  p(i,k+l)  +  p(i-l,k+l)  + 

p(i-l,k)  )  *  (  u(i,k)  +  u(i,k+l)  )  *  (  SD(i,k+l)  + 
SD (i-1 , k  +  1 )  )  -  (  p(i,k)  +  p  (i-1 ,k)  +  p (i-1 ,k-l)  + 

p (i ,k-l)  *  (  u(i,k)  +  u ( i , k-1 )  )  *  (  SD(i,k)  + 

SD ( i-1 , k )  )  )  /  (  16  *  DS  )  ;  (39) 


1  *  ", 
.■>vyy>.y»  y.  *.  v 


.2  6H 


-g  pu"  gg  (i,k)  =  {  (  p  (i,k)  +  p  (i-l,k)  * 

(  u (i,k) ) 2  *  (  HP ( i)  -  HP(i-l)  )  /  (  H(i)  *  DX)  ; 


5^  puv ( i ,  k )  =  (  (  p(i,k)  +  p(i+l,k)  )  *  u(i+l,k) 

(  v(i,k)  +  v(i+l,k)  )  -  (  p(i,k)  +  p(i-l,k)  )  * 

u(i,k)  *  (  v(i,k)  +  v ( i-1 ,k)  )  /  (  4  *  DX  ) 


^  pvs  (i,k)  =  (  (  p(i,k)  +  p(i,k+l)  )  *  SD(i,k+l) 
(  v ( i , k)  +  v ( i ,k+l )  )  -  (  p(i,k)  +  p(i,k-l)  ) 

SD(i,k)  *  (  v ( i , k )  +  v (i ,k-l)  )  /  (  4  *  DS  )  ; 


puv  ^  ( i  ,  k )  =  p  ( i , k)  *  (  u ( i , k)  +  u(i  +  l,k)  ) 


v(i,k) 


H  ( i  +  1 )  -  H(i)  )  /  (  2  *  H  *  DX  )  . 


The  remaining  variables  to  be  forecast  are  pressure 
and  density.  Density  at  the  new  time  level  is  computed 
after  (28)  is  solved  by  the  following: 


pn+1  (i,k)  =  pn  (i,k)  +  DT  *  (i, 


k)  . 


In  (44),  the  superscript  n+1  refers  to  the  forecast 
time  level  and  the  superscript  n  refers  to  the  present 
time  level.  After  the  pressure  tendency  at  the  top  of  the 
model  is  calculated  from  Eq .  (26),  the  pressure  at  the  top 

can  be  computed  by  the  following: 


p (i,kel) 


p ( i , kel )  +  DT  *  (§*) 1 


(45) 


In  the  above  finite  difference  formulas  the  vertical 
index  kel  refers  to  the  top  of  the  model.  At  this  point 
it  is  possible  to  compute  the  pressure  at  remaining  points 
using  the  hydrostatic  equation  in  the  following  finite 
difference  form: 

P (irk)  =  P(i,k+l)  +  g  *  DS  *  (HP ( i ) )  *  p(i,k)  . 

(46) 

In  the  above  finite  difference  formulas,  H  refers  to 
the  thickness  of  the  model  at  columns  where  u  is  defined, 
and  HP  refers  to  the  model  thickness  at  columns  where 
pressure,  density,  v  and  T  and  defined,  as  shown  in  Fig. 

2. 

Because  the  equation  set  used  is  hydrostatic  and 
ignores  any  buoyancy  contribution  to  vertical  motion,  a 
dry  convective  adjustment  is  made  similar  to  that 
described  by  Haltiner  and  Williams  (1980).  At  each  time 
step  the  vertical  profile  of  temperature  is  examined  for 
superadiabatic  lapse  rates  in  each  layer.  If  the  lapse 
rate  is  found  to  be  superadiabatic  in  any  layer,  it  is 
adjusted  to  a  slightly  subadiabatic  lapse  rate  by 
correcting  the  temperature  at  the  top  and  bottom  of  the 
layer.  This  process  simulates  subgrid  scale  dry 
convection  in  which  convection  transports  heat  vertically 


v 


until  a  neutral  lapse  rate  is  established.  The  kinetic 
energy  of  the  convective  transport  is  assumed  to  be 
dissipated  into  heat  energy.  The  correction  process  is 
iterative  in  that  after  a  correction  is  made  at  one  layer, 
the  adjacent  layers  are  reexamined.  The  process  is 
repeated  until  no  segment  of  the  temperature  profile  is 
superadiabatic .  Because  this  correction  is  made  at  each 
time  step,  the  resultant  temperature  change  in  each  layer 
is  small.  The  formula  used  for  the  dry-convective 
adjustment  is  as  follows: 

T1  =  T1  +  bt  ;  T0'  =  T°  -  At  ,  where 

At  =  (  (T°  -  T1)  -  0.0096  *  DS  *  H  )  /  2  .  (47) 

and  T1  and  T°  are  the  temperature  at  the  top  and  the 
bottom  of  the  layer  to  be  adjusted,  respectively.  This 
equation  approximates  the  potential  energy  conserving  form 
given  by  Haltiner  and  Williams  (1980),  and  requires  much 
less  calculation  time.  Eq.  (47)  differs  from  the  form 
given  by  Haltiner  and  Williams  (1980)  in  that  the  lapse 
rate  is  slightly  subadiabatic  after  the  adjustment  rather 
than  adiabatic.  Because  the  model  is  hydrostatic,  total 
potential  energy  is  conserved  during  the  convective 
adjustment  process  if  the  mean  temperature  of  each  layer 
is  conserved  during  the  adjustment  (Haltiner  and  Williams, 


The  x  component  of  the  wind  is  smoothed  horizontally 
at  each  time  step  using  a  simple  three  point  operator  as 
follows : 

u ( i )  =  0.90  u (i)  +  0.05  (  u(i-l)  +  u(i+l}  )  .  (48) 

This  horizontal  smoothing  is  used  to  prevent  spurious 
growth  of  short  waves  in  the  model.  The  short  waves  which 
result  in  the  model  simulations  when  smoothing  by  (48)  is 
not  performed  result  from  grid  noise  due  to  the  spatial 
averaging  in  the  finite-difference  equations.  These  short 
waves  are  minor  and  do  not  contribute  to  instability  in 
the  model,  however,  they  are  unrealistic  meteorologically 
and  are  removed  from  the  final  result  for  esthetic 


reasons . 


4. 


MODEL  VALIDATION 


Before  a  numerical  model  can  be  used,  it  is  necessary 
to  perform  basic  validation  tests  to  establish  the 
credibility  of  the  model  and  insure  that  no  errors  have 
been  made  in  programming  the  model  equations.  These  tests 
provide  an  important  means  of  debugging  the  numerical 
scheme.  The  validation  tests  performed  on  the  mesoscale 
model  consisted  of  a  simple  static  test,  calculation  of 
the  mass  continuity  and  the  kinetic  energy  budget,  and 
performing  non-linear  simulations  which  yield  physically 
expected  results  which  can  be  compared  with  the  results  of 
previous  models. 

a.  Static  test 

The  most  basic  type  of  test  which  a  numerical  scheme 
must  pass  is  a  static  test,  in  which  an  equilibrium 
initial  condition  is  allowed  to  remain  at  rest  without  any 
forcing  applied.  This  test  checks  that  all  of  the  terms 
which  should  be  zero  in  the  numerical  scheme  are  coded 
correctly  and  are  in  fact  making  no  contribution  to  the 
result.  Such  a  static  test  was  performed  with  a  domain 
size  of  30  horizontal  grid  points  and  13  vertical  grid 
points,  with  a  horizontal  grid  interval  of  45  km  and  a 
vertical  grid  spacing  of  250  m.  The  initial  state  was 
constructed  so  that  the  horizontal  derivatives  of  all 


variables  were  nonexistent  initially,  with  initial  wind  of 
zero  at  all  grid  points.  The  initial  pressure  was 
calculated  hydrostatically  by  specifying  the  surface 
temperature,  the  temperature  lapse  rate,  and  the  pressure 
at  the  top  of  the  model.  Once  the  pressure  is  calculated 
at  all  points  the  density  is  computed  using  the  equation 
of  state.  The  initial  values  of  all  variables  for  the 
static  test  are  presented  in  Table  1.  The  temperature 
lapse  rate  was  specified  at  8  K  km  As  can  be  noted, 

this  initial  distribution  of  variables  assumes  a  flat 
model,  with  no  terrain.  From  this  initial  starting 
condition  the  model  was  run  for  36  h  with  no  change  in  the 
initial  state.  An  additional  test  run  was  made  in  which 
the  bottom  of  the  model  sloped  so  that  the  value  of  H  at 
the  left  boundary  was  1  km,  while  the  top  remained  at  3 
km.  The  initial  variables  were  determined  in  the  same 
manner  as  the  previous  test,  with  the  surface  temperature 
adjusted  to  fit  the  lapse  rate  specified  as  the  surface 
sloped  upward  in  the  model.  Variables  at  the  right 
boundary  were  the  same  as  given  in  Table  1.  Because  the  s 
surfaces  in  this  case  sloped,  the  horizontal  (s) 
derivatives  no  longer  vanished,  except  at  the  top  of  the 
model.  All  pressure  surfaces  were  horizontal  with  respect 
to  sea  level.  This  test  also  showed  no  change  in  the 
initial  state  due  to  the  inclusion  of  the  terms  involving 
the  geometry  of  the  coordinate  system. 


For  all  further 


defined.  The  calculation  was  repeated  at  the  end  of  the 
simulation.  This  procedure  was  carried  out  during  the 
kinecic  energy  budget  calculations  to  be  discussed  in  the 
next  section.  The  result  after  36  h  showed  a  net  increase 
of  the  total  density  of  0.0004%.  This  small  increase  is 
negligible  and  can  be  attributed  to  truncation  errors  in 
the  finite-difference  approximations  and  computer  round 
off  error.  The  accuracy  of  the  computer  on  which  the 
model  was  run  is  limited  to  six  digits  in  single 
precision,  which  is  sufficient  to  explain  this  small 
density  variation.  These  results  indicate  that  the  mass 
continuity  of  the  model  is  satisfied,  with  no  spurious 
destruction  or  generation  of  mass  occurring  due  to 
programming  error  or  the  finite-difference  equations  used. 


c.  Testing  of  the  kinetic  energy  budget 


This  test  followed  the  general  procedure  of  the 
kinetic  energy  testing  described  by  Pielke  (1981).  Usir 
the  flux-divergence  form  of  the  equations  of  motion,  (12 
and  (13),  together  with  the  horizontal  components  of  (2) 
the  following  kinetic  energy  equation  can  be  derived. 


dt  K  = 


u  X2 
ox 
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In  this  equation  K  is  equal  to  p(u  +v  )/2.  It 
should  be  noted  that  this  kinetic  energy  equation  neglects 
the  contribution  of  s  to  the  kinetic  energy  because  it  is 
negligible  compared  to  the  horizontal  motions.  The 
equation  is  an  expression  of  the  kinetic  energy  change  per 
unit  volume  since  density  is  used  in  the  expression  for 
kinetic  energy  rather  than  mass.  By  calculating  the 
individual  terms  of  this  equation  and  integrating  over  the 
domain  of  the  model,  an  estimate  of  the  expected  kinetic 
energy  change  for  the  entire  model  can  be  made.  The 
individual  terms  in  (49)  are  calculated  using 
finite-difference  formulas  which  are  similar  to  those  used 
in  the  model  forecast.  The  finite-difference  formulas 
used  in  this  calculation  are  centered  on  grid  points  where 
the  x  component  of  the  wind  is  defined.  Performing  this 
calculation  at  each  time  step  allows  an  estimate  of  the 
rate  of  change  of  total  kinetic  energy  per  unit  volume  for 
the  model  from  which  the  total  kinetic  energy  per  unit 
volume  at  the  next  time  step  can  be  obtained.  This  is 
compared  periodically  with  the  observed  kinetic  energy, 
which  is  calculated  from  the  observed  winds  and  density, 
also  centered  on  grid  points  where  u  is  defined  so  that 
the  two  quantities  can  be  compared.  According  to  Pielke 
(1981),  if  the  two  results  closely  agree,  the  modeler  can 
be  certain  that  mistakes,  such  as  coding  errors,  are  not 
causing  significant  sources  of  unexplained  changes  of 


kinetic  energy.  To  perform  this  calculation  a  model 
domain  of  40  horizontal  grid  points  and  13  vertical  grid 
points  was  forced  by  slowly  heating  the  left  hand  side  of 
the  model  at  two  grid  columns.  The  heating  rate  was 
sufficient  to  cause  a  temperature  change  of  0.35  K  h-1  at 
the  left  hand  boundary. 

The  results  of  this  test  are  given  in  Table  2.  They 
indicate  that  the  rate  of  kinetic  energy  increase 
calculated  by  the  model  is  about  1%  h  1  less  than  that 
calculated  by  the  individual  terms  on  the  right  hand  side 
of  (46).  According  to  Anthes  and  Warner  (1978),  the 
difference  between  the  two  different  calculations  is 
partly  due  to  truncation  errors  in  the  horizontal  and 
vertical  flux  terms  in  the  equation  of  motion.  The 
remainder  of  the  difference  can  be  attributed  to  numerical 
diffusion  which  results  from  the  finite-difference  scheme 
used  in  the  forecast.  Based  on  these  results,  the 
conclusion  is  made  that  no  programming  errors  are 
contributing  to  spurious  generation  of  kinetic  energy  in 
the  model  domain,  although  some  damping  of  kinetic  energy 
does  occur.  The  initial  conditions  for  this  test  were  the 
same  as  those  given  in  Table  1.,  except  for  the  warming  of 
the  left  side  of  the  model. 


Table  2.  Results  of  kinetic  energy  budget 
calculations.  Column  A  gives  the  total  integrated  kinetic 
energy  per  unit  volume  as  observed  from  the  model 
forecast.  Column  B  gives  the  total  calculated  from 
individual  terms  in  Eg.  49. 


Hours 

A  (  J  m"3) 

B  (  J  m~3) 

1 

0.845 

0.831 

2 

2.94 

2.96 

3 

5.33 

5.43 

4 

8.00 

8.23 

5 

11.00 

11.40 

6 

14.2 

14.8 

7 

17.6 

18.6 

8 

21.2 

22.6 

9 

24.9 

26.8 

10 

28.6 

31.1 

11 

32.4 

35.5 

12 

36.1 

40.0 

13 

44.8 

14 

43.9 

49.8 

15 

48.0 

55.0 

16 

52.4 

60.7 

17 

57.1 

66.7 

d.  Nonlinear  simulations 


In  addition  to  the  basic  tests  described  above,  an 
attempt  was  made  to  duplicate  the  model  results  obtained 
by  Shen  (1980)  in  his  simulation  of  the  development  of  a 
low-level  jet.  Successful  comparison  of  such  a  simulation 
with  the  work  of  Shen  (1980)  will  give  further  support  to 
the  validity  of  this  numerical  scheme.  In  order  to 
successfully  carry  out  this  simulation,  changes  in  the 
boundary  conditions  of  the  model  were  made  as  described 


below. 


Because  Shen's  model  did  not  include  any  terrain, 


the  simulation  here  was  performed  with  a  flat  model  so 
that  the  results  could  be  compared  more  directly. 


1)  Boundary  conditions 


In  order  to  duplicate  the  work  by  Shen  (1980),  it 
was  necessary  to  change  the  boundary  conditions  from  those 
used  in  the  first  portion  of  numerical  testing.  While  the 
boundary  conditions  used  above  have  certain  advantages  for 
numerical  testing  of  the  interior  finite-difference 
formulas,  they  are  inappropriate  for  a  realistic 
simulation  of  a  low-level  jet  similar  to  that  carried  out 
by  Shen  (1980) . 

For  this  purpose  the  top  boundary  condition  was 
chosen  to  be  ^  =  0  and  ^  *  0 .  These  are  identical  to 
the  boundary  conditions  used  by  Shen  (1980) .  The  pressure 
tendency  at  the  top  of  the  model  was  set  equal  to  zero  in 
order  to  isolate  the  problem  from  the  effects  of  levels 
higher  than  3  km.  The  flux  through  the  top  of  the  model 
allows  for  the  expected  subsidence  or  downward  vertical 
motion  to  occur  when  the  model  is  forced  at  the  left 
boundary.  The  top  boundary  condition  used  in  the  earlier 
part  of  the  testing,  s  =  0,  was  also  tried  for  this 
portion  of  the  simulation,  but  caused  significant  return 
flow  at  the  upper  levels  on  the  left  hand  boundary,  which 
is  inconsistent  with  the  forcing  applied  at  this  boundary 
and  thus  considered  somewhat  unrealistic.  This  boundary 


condition  is  best  used  when  the  top  boundary  is  far 
removed  from  the  computational  region  of  interest,  which 
is  not  the  case  for  this  simulation.  For  this  reason,  and 
in  order  to  provide  better  comparison  with  the  work  of 
Shen  (1980),  the  top  boundary  conditions  were  changed  as 
described.  In  addition,  the  left  boundary  was  moved  to  a 
column  where  p  was  defined,  rather  than  a  column  of  u  as 
shown  in  Fig.  2.  The  x  and  y  components  of  the  wind  were 
calculated  next  to  the  left  boundary  by  a  simple  forward 
differencing  of  (12)  and  (13)  because  the 
finite-difference  scheme  used  for  the  interior  points 
could  not  be  used  here.  At  the  left  boundary  the  boundary 
condition  applied  was 

=  -4.44  x  10~5  (1  -  ~  )  Km"1  ,  (50) 

where  k  is  the  vertical  index.  This  provides  a 
temperature  difference  at  the  left  hand  boundary  of  2  K 
across  the  grid  interval  nearest  to  the  left  boundary. 

This  temperature  gradient  decreases  to  zero  at  the  top  of 
the  model.  This  is  consistent  with  an  assumption  of 
warming  occurring  outside  the  model  domain  which  causes 
the  left  boundary  to  be  slightly  warmer  than  the  nearest 
grid  point  in  the  interior.  This  warming  at  the  boundary 
was  also  specified  in  the  interior  portions  of  the  model 
in  a  manner  which  decreased  it  exponentially  to  zero 


farther  from  the  left  boundary  in  the  positive  x 
direction.  This  was  accomplished  by  assuming  an 
exponential  profile  of  Q  in  the  x  direction,  with  a 
maximum  value  at  the  left  end  of  the  model,  decreasing  to 
zero  in  the  x  direction  in  the  interior.  This  heating 
rate  was  a  maximum  at  the  ground  and  decreased  linearly  to 
zero  at  the  top  of  the  model.  The  surface  heating  rate  is 
given  in  Table  3.  Heating  the  model  in  this  manner 
combined  with  the  left  boundary  condition  described  above 

Table  3.  Surface  heating  rate  (  K  s  ^  x  10  )  . 

Values  decrease  linearly  to  zero  at  the  top  of  the  model. 


A  X 


Q/C 


forcing  mechanisms  is  similar  to  that  used  by  Shen  (1980)  , 
who  specified  a  heating  rate  at  the  left  hand  boundary 
only.  For  this  reason  Shen ' s  model  became  considerably 
warmer  at  the  left  hand  boundary  than  at  the  nearest 
internal  grid  point,  building  a  horizontal  temperature 
gradient  of  20  K  over  the  first  grid  spacing  adjacent  to 
the  left  boundary  over  his  28  h  forecast  period.  The 
forcing  used  here  implies  a  strong  horizontal  diffusion  of 
heat,  and  prevents  a  large  temperature  gradient  from 
occurring  at  the  left  boundary  as  in  Shen ' s  model.  The 
pressure  gradient  at  the  right  boundary  was  assumed  zero, 
following  the  procedure  of  Shen  (1980).  The  above  forcing 
mechanism  was  chosen  because  it  is  similar  to  that  used  by 
Shen  (1980),  although  it  is  not  physically  realistic 
because  no  diurnal  variation  of  the  heating  is  allowed. 

2)  Rederivation  of  Richardson  1 s  equation . 

Because  the  form  of  Richardson's  equation  used  in  the 
earlier  validation  testing  was  derived  using  a  top 
boundary  condition  of  s  =  0 ,  it  was  necessary  to  rederive 
this  equation.  The  following  equation  used  in  the 
derivation  is  presented  by  Djuric  and  Das  (1982)  : 

ds  1__  d£  Q  _1 

=  *  |p  dt  +  C  T  ‘  h  7  '  •  (51) 

ds 

Applying  the  top  boundary  condition  ^  =  0  and 


assuming  Q  =  0  at  the  top  of  the  model  gives  the  following 


relationship : 


(£Ej 
dt  1 


=  -  (  XE  v  •  HV)  1 


Also  given  by  Djuric  and  Das  (1982)  is  the  following 


equation : 


=  (v  •  VP) 1  -  g  f  v  •  HpV  ds '  +  (|f ) ±  -  gp1Hs1  , 


which  can  be  written  as 


<§?>!  =  (V  •vpl1  -  gp^s,  ,  (If), 


since  the  integral  vanishes  at  s=l. 

Combining  (52)  and  (54)  gives  the  following: 


(|f)x  =  -  *gl  (V  •HV)1  -  (V  •Vp)1  +  gp1Hs1 


Substituting  (53)  into  (51)  yields: 


§f  =  ~  ~  V  •  Vp  -  g  (V  ■  HpV) ds 1  -  gp1Hs1 


+  -2_  -  i  v  .  h\^ 

CpT  H 


Eq.  (55)  can  be  substituted  into  this  equation  to 


6s _ 1 

6s  TP 


V  •  Vp  - 


-  Xgl  (V  •  HV)  x 


(V  *  Vp)  1  -  g 


1 

(V  •  HpV) ds ’ 


+ 


Q 

C  T 
P 


s 


1 

H 


V  •  HV 


(57) 


which  can  be  integrated  to  yield  the  Richardson's  equation 
as  follows: 


1  -f  'V.vpij  , 


0 


D  -  $  •  Vp  Q 

TP  +  cpT 


£  V  -  HV)  ds  ' 
ri 


o 


+  I1  <v  hv^i  /  H 

0 


d_s  ' 

P 


:58) 


where  the  function  D  is  the  same  as  previously  defined. 
Application  of  the  other  top  boundary  condition  eliminates 
the  first  term  in  this  equation.  It  can  be  noted  that 
only  one  term  of  this  form  of  Richardson's  equation  is 
different  from  the  previous  form.  For  this  reason  little 
additional  programming  was  necessary  to  include  the  new 
top  boundary  conditions.  The  equation  for  density 
tendency  becomes 


V  •  V  p 


+  -fi 
TP 


6t  1 


ghp1s1 


P 


D  -  V  •  V  p  + 

TP 


v 


(59) 


Derivation  of  this  equation  follows  that  given  in  the 


Appendix,  with  the  term  involving  s  at  the  top  included. 


It  can  be  noted  from  (55)  that  1=  ®  requires  the 

following  relationship: 

(*§)j_  (V  •  HV)X  =  gp1Hs1  _  (60) 

This  equation  is  used  to  determine  u  at  the  top 
boundary  by  integrating  from  the  right  boundary  where  u  is 
zero,  to  the  left  boundary.  This  provides  a  slight 
horizontal  convergence  at  the  top  of  the  model  to 
counteract  the  subsidence  so  that  the  resultant  pressure 
tendency  at  the  top  of  the  model  vanishes.  This  is 
necessary  for  the  wind  field  at  the  top  of  the  model  to  be 
dynamically  consistent  with  the  boundary  conditions 
stated.  Using  (60)  to  find  u  at  the  top  boundary  results 
in  a  very  small  u  field  because  the  necessary  convergence 
required  by  (60)  is  small.  Following  the  integration  of 
(60),  the  v  field  at  the  top  is  obtained  using  the 
horizontal  equation  of  motion  for  v,  with  the  assumption 
that  the  vertical  advection  term  is  small.  The  resultant 
correction  of  u  and  v  at  the  top  of  the  model,  while 
dynamically  consistent  with  the  boundary  conditions, 
actually  has  very  little  effect  on  the  results  obtained. 
This  was  shown  by  making  an  additional  simulation  in  which 
the  simple  condition  u  =  v  =  0  at  the  top  was  used, 
effectively  ignoring  the  right  hand  of  (60).  While  not 
dynamically  consistent,  virtually  identical  results  were 
obtained.  This  is  due  to  the  fact  that  the  pressure 


tendency  at  the  top  is  quite  important,  appearing  in  the 
density  tendency  equation  at  all  levels,  while  the  wind  at 
the  top  of  the  model  is  used  only  in  the  vertical 
advection  terms  in  the  equations  of  motion  and  in  the 
calculation  of  internal  friction  at  the  top  of  the  model 
where  a  vertical  derivative  of  the  wind  is  required.  Both 
the  vertical  advection  terms  and  the  internal  friction 
terms  involved  are  very  small,  and  have  little  effect  on 
the  rest  of  the  computation. 

3 )  Initial  conditions 

The  initial  conditions  used  for  the  simulations 
consisted  of  a  slightly  stable  boundary  layer,  capped  by 
an  inversion.  Initial  values  for  all  variables  are  given 
in  Table  4.  The  temperature  lapse  rate  assumed  in  the 
boundary  layer  was  8  K  km-1,  slightly  less  than  dry 
adiabatic.  The  pressure  surfaces  were  all  calculated 
using  the  hydrostatic  equation  based  on  the  initial 
temperature  profile  and  the  pressure  at  the  top  of  the 
model.  Density  was  calculated  from  the  pressure  and 
temperature  using  the  equation  of  state.  The  initial  wind 
field  was  zero  throughout.  The  initial  values  of  all 
variables  were  horizontally  uniform,  because  the  model  was 
assumed  to  be  flat. 


Results  of  the  simulation 


The  equations  were  integrated  for  28  h  beginning  with 
the  initial  conditions  given  in  Table  4.  The  resulting 
flow  was  forced  by  a  combination  of  the  left  boundary 
condition  and  the  heating  in  the  interior  of  the  model. 
Results  are  shown  in  the  sets  of  Figs.  3-9.  The  first 
figure  in  each  set  consists  of  an  x  -  s  cross  section 
showing  an  isotach  analysis  of  the  v  component  of  the 

Table  4.  Initial  values  of  variables  for  nonlinear 
simulation  of  low-level  jet. 


Height  (m) 


P  (mb)  Density  (kg/m  ) 


267.0 

677.0 

0.898 

269.0 

699.0 

0.920 

271.0 

721.5 

0.942 

273.0 

744.6 

0.965 

275.0 

768.3 

0.989 

277.0 

792.5 

1.012 

264.0 

817.3 

1.096 

266.0 

844.1 

1.124 

268.0 

871.7 

1.151 

270.0 

899.9 

1.180 

272.0 

928.8 

1.209 

274.0 

958.4 

1.238 

988.7 
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Fig.  4a.  Isotach  analysis  at  8  h  for  v. 
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Fig.  4b.  Isotherms  at  8  h. 
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Fig.  5a.  Isotach  analysis  at  12  h  for  v. 
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Fig.  5b.  Isotherms  at  12  h. 
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Fig.  8a.  Isotach  analysis  at  24  h  for  v. 
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Fig.  8b.  Isotherms  at  24  h. 


wind,  while  the  second  figure  in  each  set  consists  of  an 
x  -  s  cross  section  showing  an  isotherm  analysis. 

Selected  isotachs  in  the  first  figure  of  each  set  are 
labelled  in  m  s  1.  Selected  isotherms  in  the  second 
figure  of  each  set  are  labelled  in  K.  The  sets  of  figures 
are  given  at  intervals  of  4  h,  beginning  with  the  results 
at  4  h. 


1 )  Propagation  of  the  low- leve 1  jet 

Because  the  coriolis  parameter  f  was  considered 
constant  during  the  calculation,  the  x  axis  will  be 
considered  an  axis  of  constant  latitude  in  the  following 
discussion.  For  this  reason,  negative  values  of  u  are 
referred  to  as  easterly,  while  positive  values  of  v  are 
referred  to  as  southerly  in  the  following  comments.  As 
shown  by  the  isotach  analysis  in  Fig.  4,-  a  low-level 
southerly  wind  develops  adjacent  to  the  left  boundary 
where  the  maximum  warming  occurs.  The  maximum  value  of 
the  wind  remains  near  the  left  boundary,  with  the 
southerly  wind  speed  increasing  throughout  the  integration 
period.  This  result  is  similar  to  that  obtained  by  Shen 
(1980),  with  the  exception  that  the  maximum  wind  value  did 
not  show  the  slight  eastward  movement  obtained  by  Shen 
(1980).  The  maximum  value  of  the  wind  obtained  is  less 
than  one-half  the  value  obtained  by  Shen  (1980).  This  is 
due  to  the  difference  in  forcing  mechanisms  between  the 


two  models  discussed  earlier.  The  temperature  gradient 
achieved  by  the  heating  mechanism  in  this  model  was  only 
about  7  K  across  the  entire  domain.  The  temperature 
increase  at  the  lowest  level  on  the  left  hand  boundary  was 
about  7  K,  less  than  the  heating  rate  obtained  by  Shen. 
Because  of  the  large  difference  in  horizontal  temperature 
gradients  obtained  in  the  two  models,  it  is  expected  that 
the  maximum  winds  obtained  would  be  less  in  the  present 
study.  While  the  maximum  values  obtained  are  lower,  the 
values  of  southerly  wind  produced  in  the  interior  of  the 
model  are  similar,  as  shown  by  Fig.  10.  Values  for  the 
maximum  southerly  wind  obtained  by  Shen  (1980)  presented 
in  Fig.  10  are  evaluated  from  the  isotach  analysis 
presented  in  his  paper.  The  numerical  scheme  used  in  the 
present  study  causes  the  pressure  fall  to  propagate 
rapidly  across  the  model  domain  when  the  model  is  heated 
at  the  boundary.  Attempts  were  made  to  achieve  horizontal 
temperature  gradients  similar  to  that  obtained  by  Shen 
(1980)  by  setting  Q  =  0  in  the  interior  of  the  model  and 
heating  strongly  at  the  boundary,  but  it  was  not  possible 
to  achieve  as  strong  a  pressure  gradient  across  a  small 
grid  interval  as  did  Shen  (1980)  without  excessive 
heating . 

2 )  Subside nee -type  inversion 


Part  b  in  Figs.  3-9  shows  the  distribution  of 


TIME  (H) 


Fig.  10.  Comparison  of  maximum  speed  of  southerly 
wind  400  km  from  left  hand  boundary  in  Shen ' s  model  with 
TAMU  two-dimensional  model. 

temperature  at  4  h  intervals,  with  isotherms  labeled  in  K. 
These  results  show  the  time  evolution  of  a  subsidence  type 
inversion  very  similar  to  that  observed  by  Shen  (1980)  in 
his  study.  Consistent  with  his  results,  the  subsidence 
inversion  layer  moves  downward  with  time,  with  the 
greatest  vertical  displacement  occurring  at  the  left 
boundary  and  lesser  vertical  displacement  occurring 
eastward  across  the  model  domain.  The  strength  of  the 
inversion  is  also  noted  to  decrease  slightly  in  the 
positive  x  direction.  As  noted  by  Shen  (1980) ,  the 
occurrence  of  a  subsidence  inversion  and  its  importance  to 
the  low-level  jet  have  been  noted  by  many  investigators, 


including  Bonner  (1966),  Beckman  (1973),  and  Damiani 
(1979) .  It  is  significant  that  the  subsidence  was 
sufficient  to  prevent  the  inversion  from  disappearing 
despite  the  differential  vertical  heating  applied.  Figure 
11  shows  the  distribution  of  s  at  28  hours  in  cm  s-1.  The 
result  of  subsidence  occurring  throughout  the  model  domain 
i"  at  first  surprising  in  light  of  the  fact  that  Q  makes  a 
positive  contribution  to  Richardson's  equation.  Scale 
analysis  of  the  individual  terms  involved  shows  that  near 
the  level  of  the  maximum  wind,  the  term  — ^  V  -  1/  H  is 
about  1  order  of  magnitude  larger  than  the  other  terms  in 
the  second  integrand.  This  term  makes  a  negative 
contribution  to  vertical  motion  due  to  the  horizontal 
divergence  of  the  wind  field.  The  last  term  is  also 
negative  due  to  the  convergence  which  occurs  at  the  top  of 
the  model  to  satisfy  the  top  boundary  condition,  and  is 
comparable  in  magnitude  to  the  other  divergent  term 
discussed  above.  The  remaining  terms  make  a  positive 
contribution  which  is  too  small  to  offset  the  contribution 
of  the  negative  terms.  Setting  u  =  0  at  the  left  boundary 
allows  the  development  of  a  positive  vertical  motion  in 
the  areas  where  heating  is  applied  because  of  the 
horizontal  convergence  which  occurs  near  the  left  boundary 
as  the  x  component  of  the  wind  is  forced  to  slow  down. 
Figure  11  indicates  that  mass  is  entering  the  top  of  the 
model  to  partially  compensate  for  the  mass  flux  through 


90 


180  270  360  450  540  630  720  810 

X(  KM) 

Fig.  11.  Vertical  motion  s  at  28  h  in  cm  s-1. 


the  left  boundary,  resulting  in  subsidence  throughout  the 
model  domain.  Figure  12  shows  the  circulation  in  the  x-s 
plane  given  by  the  stream  function  defined  by 


u  * 


(61 


s 


1_ 

H  <bx 


(62 


The  stream  function  values  in  Fig.  12  were 
approximated  by  integrating  (61)  from  the  surface,  where 
the  stream  function  is  arbitrarily  equal  to  zero,  to  the 
top  of  the  model.  A  more  accurate  stream  function 
calculation  would  involve  solving 

,  (63 
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Fig.  12.  Streamf unction  in  m  s 
where  f  .«$*  - -I  §|  •  ,64, 

Solution  of  (63)  would  require  boundary  conditions 
for  on  all  boundaries.  The  approximate  solution  given 
by  (61)  neglects  the  nondivergent  part  of  the  circulation 
in  the  vertical  plane.  Figure  12  also  shows  that  mass  is 
entering  the  top  of  the  model  and  exiting  the  lower  left 
boundary,  in  agreement  with  Fig.  11.  Values  of  s 
estimated  from  Fig.  12  by  applying  (62)  are  in  fair 
agreement  with  the  calculated  values,  giving  some  further 
justif ication  for  the  integration  of  (61)  rather  than  the 
complete  solution  of  (63). 


3)  The  relation  between  the  low-level  jet  and  the 
inversion. 

As  shown  by  comparison  of  the  isotach  and  isotherm 
analyses  in  Figs.  3-9,  the  maximum  value  of  the  southerly 
wind  occurs  just  at  or  below  the  inversion  base  in  each 
frame.  By  28  h  the  bottom  of  the  inversion  layer  is 
almost  coincident  with  the  maximum  southerly  wind  band. 

As  shown  in  the  isotach  analyses,  the  level  of  the  maximum 
southerly  wind  increases  eastward  slightly  as  the  height 
of  the  inversion  increases  eastward  spatially.  This 
result  also  was  obtained  by  Shen  (1980)  .  For  comparison, 
a  model  simulation  was  also  carried  out  in  which  the 
initial  conditions  were  similar  with  the  exception  that 
the  initial  temperature  profile  was  slightly  stable 
throughout  with  no  inversion  layer.  The  isotachs  in  the 
x  -  s  plane  for  the  southerly  component  of  the  wind 
achieved  at  28  h  are  presented  in  Fig.  13  in  m  s  1. 
Comparison  of  these  isotachs  with  those  in  Fig.  9  b  shows 
that  the  resulting  wind  is  deeper  in  the  no-inversion 
case,  with  a  weaker  vertical  wind  shear  produced  above  the 
level  of  the  maximum  wind.  The  strength  of  the  maximum 
wind  obtained  is  similar  due  to  the  similar  forcing 
mechanism.  The  only  difference  which  reasonably  could  be 
expected  concerns  the  vertical  variation  of  the  wind, 


Because  pressure  is  not  allowed  to  vary  along  the  y 
axis  in  this  study,  the  u  component  of  the  wind  is  the 
isallobaric  wind.  Figure  14  shows  the  time  evolution  of 


the  maximum  values  of  the  easterly  and  southerly  wind. 

These  results  shov/  that  the  easterly  wind  is  stronger 
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Fig.  14.  Time  variation  of  the  maximum  southerly  and 
easterly  wind. 

for  the  first  6  h  of  the  simulation,  after  which  the 
southerly  wind  becomes  increasingly  stronger  with  time  due 
to  the  Coriolis  effect.  The  easterly  wind  appears  to 
become  relatively  steady  after  16  h.  These  results 
indicate  that  the  southerly  wind  becomes  increasingly 
quasi-geostrophic  with  time  during  the  simulation.  The 
easterly  wind  does  not  become  increasingly  quasi- 
geostrophic  because  this  would  require  the  development  of 
a  pressure  gradient  in  the  y  direction  which  is  not 
allowed  in  this  two-dimensional  simulation.  Examination 
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of  Figs.  3-9  shows  that  the  thermal  wind  relationship  is 
approximated  above  the  level  of  the  maximum  wind,  above 
which  the  wind  decreases  with  height.  This  would  result 
in  a  thermal  wind  at  these  levels  whose  magnitude  is  less 
than  zero.  This  is  consistent  with  the  warming  at  the 
left  boundary.  '  Below  the  level  of  the  maximum  wind  the 
surface  frictional  effects  do  not  allow  the  thermal  wind 
to  be  approximated.  These  results  are  consistent  with  the 
gecstrophic  adjustment  reported  by  Shen  (1980)  in  his 
simulation . 


5. 


CONCLUSIONS  AND  RECOMMENDATIONS 


As  a  result  of  this  study,  the  following  conclusions 
can  be  made : 

1.  Conservation  of  mass  is  not  violated  by  the 
finite-difference  equations  used  at  interior  grid  points 
in  this  model.  Any  loss  or  gain  of  mass  which  occurs  in 
the  model  domain  is  due  to  differential  flux  of  mass 
through  the  boundaries  of  the  model. 

2.  Kinetic  energy  production  by  the  model  is 
slightly  lower  than  that  predicted  by  an  analysis  of 
individual  terms  in  the  kinetic  energy  budget.  This 
analysis  shows  that  no  spurious  kinetic  energy  is  created 
by  the  finite-difference  equations. 

3.  The  finite-difference  scheme  is  numerically 
stable  for  at  least  28  h  with  the  parameters  used  in  this 
study. 

4.  The  results  obtained  in  a  simulation  of  the 
development  of  a  low-level  jet  are  very  similar  to  those 
obtained  by  Shen  (1980).  The  low-level  jet  showed 
comparable  development  and  was  associated  with  an 
inversion  which  capped  the  maximum  southerly  wind.  The 
inversion  was  maintained  by  subsidence  in  the  model.  The 
quasi-geostrophic  adjustment  of  the  southerly  flow  was 
similar  in  both  cases.  The  model  results  differed  in  the 
strength  of  the  maximum  wind  due  to  the  different  pressure 


gradients  which  were  produced  as  a  result  of  differential 
heating.  A  control  case  involving  no  inversion  showed  a 
deeper  wind  field  with  weaker  vertical  wind  shear  above 
the  level  of  maximum  wind. 

5.  The  finite-difference  equations  used  in  this 
study  were  similar  to  those  used  in  the  three-dimensional 
model  of  Djuric  and  Das  (1982),  except  for  those  used  in 
the  advection  scheme.  In  light  of  the  fact  that  no 
numerical  instabilities  occurred  in  the  two-dimensional 
version  used  here,  it  is  likely  that  the  difficulties 
experienced  in  the  three-dimensional  model  of  Djuric  and 
Das  (1982)  were  due  to  unresolved  programming  errors, 
possibly  related  to  the  three-dimensional  advection 
scheme.  The  three-dimensional  model  used  an  enstrophy  and 
kinetic  energy  conserving  advection  scheme  (Djuric  and 
Das,  1982),  which  was  not  used  in  the  simplified  model 
presented  here  because  it  is  not  appropriate  for  a 
two-dimensional  model.  The  finite-difference  forms  of 
Richardson's  equation  and  the  density  tendency  equation 
used  here  appear  to  be  sound  because  they  were  used 
successfully.  A  three-dimensional  model  incorporating 
forms  of  Richardson’s  equation  and  the  density  tendency 
equation  similar  to  those  used  here  appears  feasible  if 
the  advection  scheme  used  is  carefully  formulated. 

The  following  recommendations  are  made: 

1.  The  model  should  be  improved  by  adding  moisture 


to  the  equations  used,  with  careful  attention  given  to  the 
convective  adjustment  process. 

2.  Improvements  should  be  made  in  the  internal 
friction  terms  to  provide  a  more  realistic  simulation  in 
the  boundary  layer.  Also,  resolution  should  be  increased 
in  the  lower  boundary  layer. 

3.  Further  testing  should  be  undertaken  by 
incorporating  real  data  into  the  model  to  attempt  to 
simulate  actual  atmospheric  processes. 

4.  Further  attempts  at  modeling  of  a  low-level  jet 
should  use  more  realistic  lateral  boundary  conditions  than 
those  used  here.  In  addition,  the  forcing  mechanism 
should  be  more  physically  realistic  than  that  used  here  or 
by  Shen  (1980).  Diurnal  variations  of  heating  should  be 
included.  Such  simulations  should  be  three-dimensional  if 
possible . 

5.  The  model  should  be  extended  to  three  dimensions, 
with  careful  attention  given  to  the  advection  scheme  used. 
The  validation  tests  used  here  should  be  repeated  on  the 
new  three  dimensional  version,  followed  by  verification 
tests  involving  real  data.  Succesful  completion  of  such  a 
model  would  fulfill  the  requirements  of  Djuric  and  Das 
(1982)  given  previously. 
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APPENDIX 


The  following  derivations  are  given  by  Djuric  and  Das 
(1982)  and  are  repeated  here  for  convenience. 

a.  Derivation  of  the  pressure  tendency  equation  (7) 

Differentiation  of  the  hydrostatic  equation  (3)  with 
respect  to  time  yields: 


5i  5t  =  "  g  H  It 


(65) 


Elimination  of  between  (65)  and  (4)  gives: 


6_  6 

bs 


5t  ’  -g  H  (  -  V  •  pv  -  ps  -  f  V  •  VH  ) 


(66) 


which  can  be  written 


SI  =  g  V  •  Hptf  ♦  gH  ps 


(67) 


Integration  of  (67)  from  some  level  s  to  the  top  of 
the  model  gives 
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6t ' s  dt  1  g 


V  •  HpV  ds  +  gHps  -  gHp1s1 


(68; 


which  is  (7) . 


The  First  Law  of  Thermodynamics  in  entropy  form  is 
Q  1_  dT  _R  d£ 

C  T  '  T  dt  '  pC  dt  .  (71 

?  P 


Writing  the  equation  of  state  (6)  as 
dT  _  dp  dp 

T  "  p  “  p  (72 

and  substituting  into  (71)  gives  the  result 

Q  1  dp  1  dp  R  dp 

C  T  p  dt  p  dt  pC  dt  (73 
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or 
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This  can  be  written  as 

_1  dp  _  1_  d£  Q 

p  dt  =  |*p  dt  C  T  ,  (75 

It 


where  y  =  — p 

Cv 

Eliminating  between  (70)  and  (75)  and  solving  for 


&s 


Bra 


Adding  V  •  Vp  +  s  to  both  sides  of  (68)  gives 

1 

|£  +  V-Vp  +  s||  =  V-Vp  +  s  If  -  g  J  V  •  HpV  ds  + 


gHps  -  gHp1s1  +  (|^)x 


The  hydrostatic  equation  can  be  used  to  eliminate  the 
bracketed  terms  since  they  cancel  one  another. 

Substituting  the  result  into  (77)  gives 

M  =  _  VP  +  .  I  (&E)  +  _0_  _  1  v  . 

ds  TP  TP  TP  dP  1  cdt  h  '  ,79) 


where  D  =  g  V  •  Hp^/  ds'  and  gp^Hs^  =  0  because  =  0. 


Integration  of  (79)  gives 


D  -  V  •  Vp  Q  1  _  4  .  ,  1  .bp.  ds' 

'  TP  CpT  H  V  *  )  ds  T  bt)  l  J  P 
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which  is  Richardson's  equation  (9) 


c.  Derivation  of  equation  for  pressure  tendencv  at  s=l 


Application  of  the  top  boundary  condition  s  =  0  at  s 


1  to  (80)  and  rearranging  terms  gives 


which  is  (10)  . 


d.  Derivation  of  density  tendency  equation . 


Eq.  (75)  can  be  written 


||,^.Vp,4|£  =  e_  |'«£*tf.Vp  +  s|§ 


-  e£_ 
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Substituting  (68)  into  (82)  and  simplifying  using 
the  hydrostatic  equation  and  the  top  boundary  condition 
gives 


It  =  -  V  •  7p  -  s  ||  +  "  P 


D  -  V  •  Vp  _Q_ 


TP 


C  T, 
P  J 

(83) 


which  is  Eq.  (8). 
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